AFML-TR-67-121 


Part  IV 

07/36  3^  OFFICIAL  FILE  COPY 


EVALUATION  OF  MOLECULAR  WEIGHT  FROM 
EQUILIBRIUM  SEDIMENTATION 

PART  IV:  APPLICATION  OF  TIKHONOV'S  REGULARIZING  FUNCTIONS 


MAT  ATI  A  HU  T.  GEHATIA 
DONALD  R.  WIFE 


TECHNICAL  REPORT  AFML-TR-67-121,  PART  IV 


AUGUST  1970 


This  document  has  been  approved  for  public  release 
and  sale;  its  distribution  is  unlimited. 


BEST  AVAILABLE  COPY 


AIR  FORCE  MATERIALS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO 

ADOWJOB*© 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used  for  any  purpose 
other  than  in  connection  with  a  definitely  related  Government  procurement  operation, 
the  United  States  Government  thereby  incurs  no  responsibility  nor  any  obligation 
whatsoever;  and  the  fact  that  the  government  may  have  formulated,  furnished,  or  in 
any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not  to  be  regarded 
by  implication  or  otherwise  as  in  any  manner  Licensing  the  holder  or  any  other  person 
or  corporation,  or  conveying  any  rights  or  permission  to  manpfacture,  use,  or  sell  any  . 
patented  Invention  that  may  in  any  way  be  related  thereto. 

This  document  has  been  approved  for  public  release  and  sale;  its  distribution 
is  uni imi ted. 


Copies  of  this  report  should  not  be  returned  unless  return  is  required  by  security 
considerations,  contractual  obligations,  or  notice  on  a  specific  document. 


iU0-O<  tub.  19T0-<.;O'.0S-7-71-  12,^ 


o  C 


AFML-TR-67-121 
Part  IV 


August  1970 


ERRATA  -  October  1970 

The  following  corrections  are  applicable  to  AFML-TR-67-121,  Part  IV, 
(UNCLASSIFIED  Title)  Evaluation  of  Molecular  Weight  From  Equilibrium 
Sedimentation.  Part  IV;  Application  of  Tikhonov’s  Regularizing  Functions. 
UNCLASSIFIED  Report,  August  1970: 


Page  6 


Change  equation  (8)  to  read 


0 


Air  Force  Materials  Laboratory 
Air  Force  Systems  Command 


Wright-Patterson  Air  Force  Base,  Ohio 


AFML-TR-67-121 
Part  IV 


EVALUATION  OF  MOLECULAR  WEIGHT  FROM 
EQUILIBRIUM  SEDIMENTATION 

PART  IV:  APPLICATION  OF  TIKHONOV'S  REGULARIZING  FUNCTIONS 


MATATIAHU  T.  GEHATIA 
DONALD  R.  WIFE 


This  document  has  been  approved  for  public  release 
and  sale;  its  distribution  is  unlimited. 


AFML-TR-67-121 
Part  IV 


FOREWORD 


This  report  was  prepared  by  the  Polymer  Branch  of  the  Nonmetal  lie 
Materials  Division.  The  work  was  initiated  under  Project  No.  Tiki, 
"Fundamental  Research  on  Macromolecular  Materials  and  Lubrication 
Phenomena,"  Task  No.  73^203>  "Fundamental  Principles  Determining  the 
Behavior  of  Macromolecules,"  with  Dr.  M.  T.  Gehatia  acting  as  task 
scientist.  The  work  was  administered  under  the  direction  of  the  Air 
Force  Materials  Laboratory,  Directorate  of  Laboratories,  Air  Force  System 
Command,  Wright-Patterson  Air  Force  Base,  Ohio  45^33. 

The  report  covers  research  conducted  from  January  I969  to 
February  1970.  The  manuscript  was  released  by  the  authors  In  April  1970 
for  publication  as  a  technical  report. 

This  technical  report  has  been  reviewed  and  is  approved. 


WILLIAM  E.  GIBBS 
Chief,  Polymer  Branch 
Nonmetal  lie  Materials  Division 
Air  Force  Materials  Laboratory 


AFML-TR-67-121 
Part  IV 


TABLE  OF  CONTENTS 

SECTION  PAGE 

I  INTRODUCTION  1 

I  I  IMPROPERLY  POSED  PROBLEM  3 

III  COMPUTATION  8 

IV  RESULTS  OF  COMPUTATION  11 

V  CONCLUSIONS  I5 

APPENDIX  COMPUTER  PROGRAM  LISTINGS  25 

REFERENCES  34 


V 


ILLUSTRATIONS 


AFML-TR-67--121 
Part  IV 


FIGURE  page 

1.  Calculation  of  the  inverse  problem,  f(S)  distribution 
without  regularization.  The  solid  curve  is  the  original 
distribution  f,Cs)»  the  dashed  curve  is  the  resulting 
distribution  (s). 

2.  Minimizing  the  smoothing  functional  [f  (s)  jti  ( ^  )  ] ,  A 

plot  of  loga  versus  M®.  17 

3.  Determination  of  MWD  for  a  mono-modal  distribution  by 

applying  Tikhonov's  regularizing  function  and  the 
concentration  kernel.  The  solid  curve  represents  the 
original  f(s)^and  the  crosses  represent  the  computed 
distribution  f°*  (s) .  ]8 

4.  Determination  of  MWD  for  a  mono-modal  distribution  by 
applying  Tikhonov's  regularizing  function  and  the 
concentration  kernel.  TJie  solid  gradient  curve 
represents  the  original  f(t)  and  the  crosses  represent 

the  computed  distribution  f°'(s).  19 

5.  The  ultracentrifugal  u(^  )  functions  computed  from  the 
original  distribution  and  the  calculated  distribution. 

The  solid  curve  represents  the  original  ir(C)‘  The 

dotted  line  is  [u‘(  ^  )-u*  (  ^  )  ]  x  103.  20 

6.  Evaluation  of  the  bi-modal  distribution.  The  solid 
line  represents  the  original  distribution.  The  dotted 
line  was  the  resulting  distribution  using  Simpson's 
formula  for  integration.  The  dashed  curve  resulted 

from  application  of  trapezoidal  integration.  21 

7.  Resulting  distributions  evaluated  according  to 
different  error  criterion.  The  solid  curve  is  the 

original  distribution.  22 

8.  Evaluation  of  a  bi-model  distribution  for  various 

values  of  (=  Xm  ).  The  original  distribution 

max 

is  Indicated  by  the  solid  curve.  23 

9.  Evaluation  of  a  b I -modal  distribution  using  a  modified 
concentration  kernel.  Here  ^  =  (r^-r,^)/ (r,  ^-r  2)^ 

3  0  3 

The  solid  curve  represents  the  original  distribution.  24 


vl 


AFML-TR-67-121 
Part  IV 

SECTION  I 

INTRODUCTION 


The  determination  of  molecular  weight  distribution  (MWD)  from 
equilibrium  sedimentation  was  investigated  in  1928  by  Rinde  (Reference  1). 
Fujita  (Reference  2)  expressed  the  relationship  between  the  frequency 
function,  f  (m)  ,  of  MWD  and  an  experimental  curve  u(^),  obtained  from 
ultracentrifugal  sedimentation  in  the  form  of  the  two  following 
equat i ons : 


mox 


Xme 


-Xm^ 


I-  e 


rXm 


f (m) dm 


u(C)  = 


I  dC(C) 
Cq  dC 


-I 


■"max 


I  -  e 


-Xm 


f  (m)  dm 


(I  ) 


(2) 


where  C  is  concentration;  Co,  initial  concentration;  m,  molecular  weight, 
r,  distance  from  the  center  of  rotation;  to  the  meniscus  and 

bottom  respectively;  p  ,  density  of  solution;  V,  partial  specific  volume 
of  polymer;  cu  ,  angular  viscosity;  R,  universal  gas  constant;  and 
T,  temperature  in  °K;  and 


X  = 


2RT 


The  determination  of  MWD  can  now  be  defined  as  a  solution  of  Equation  1 
or  2. 
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In  recent  years  the  need  to  know  the  MWD  of  polymers  has  significantly 
Increased  and  spurred  many  new  attempts  to  develop  a  good  computational 
technique  for  this  purpose  (References  3  through  17).  The  determination 
of  MWD  was  also  an  objective  of  this  laboratory  (References  18  through  21). 
Unfortunately,  the  efforts  to  develop  a  good  method  for  MWD  determination 
were,  in  general,  not  very  successful  and  the  results  proved  to  be 
disappointing.  In  some  cases  such  results  lacked  precision  and 
reproducibility  and  they  did  not  seem  to  provide  a  routine  which  would 
work  in  all  circumstances.  Elucidation  of  this  enigmatic  behavior  of 
FuJIta's  equations,  and  the  difficulties  arising  during  their  solution 
is  a  subject  of  Investigation  in  this  laboratory. 

Lee  (Reference  21)  first  recognized  that  Equations  1  and  2  constitute 
an  Improperly  Posed  Problem  (IPP)  in  the  Hadamard  sense  (Reference  22) 
and,  therefore,  a  conventional  mathematical  approach  cannot  provide 
satisfactory  solutions.  There  is,  however,  the  possibility  of  finding 
an  approximate  solution  to  IPP  in  general  and  to  FuJIta's  equations  in 
particular  by  applying  Tikhonov's  regularizing  functions.  Application 
of  this  method  and  solution  of  Fujita's  equations  is  the  objective  of 

this  Technical  Report. 

\ 
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SECTION  II 

IMPROPERLY  POSED  PROBLEM 

The  notion  of  a  mathematical  IPP  (Improperly  Posed  Problem,  also 
known  as  an  111  Posed  Problem  or  Incorrectly  Formulated  Problem)  goes 
back  to  Hadamard  (Reference  22)  in  conjunction  with  the  Cauchy  problem 
for  the  Laplace  equation.  At  that  time  Hadamard  believed  an  IPP  Is 
evidence  that  the  expression  under  consideration  does  not  have  any  real 
physical  meaning,  and  therefore,  it  cannot  express  any  physical  rule  or 
situation.  This  belief  proved  to  be  wrong,  and  the  existence  of  IPP  was 
established  in  many  other  "physical"  expressions  such  as  the  heat 
equation  for  negative  time  and  Cauchy  data  on  the  boundary,  the  non- 
hyperbolic  Cauchy  problem  for  the  wave  equation,  the  Inverse  problem  of 
potential,  and  a  number  of  inverse  problems  for  differential  and  Integral 
equations. 

In  the  recent  decade  IPP  was  intensively  investigated.  The  following 
considerations  are  based  upon  the  works  of  C.  Pucci  (References  23 
through  28);  F.  John  (References  29  through  32);  M.  M.  Lavrentiev 
(References  33  through  43);  A.N.  Tikhonov  (References  44  through  61); 
V.K.  Ivanov  (References  62  through  66);  V.A.  Morozov  (References  67 
through  70) ;  and  others  (References  71  through  83) • 

Consider  the  equation: 

u(^)  =  Af(m)  =  Q[C,f(m)| 
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where  A  is  an  operator,  ueU,  fc  F,  U,  and  F  are  some  complete  metric 
spaces,  and  Af  is  a  function  with  the  domain  of  definition  F  and  the 
range  of  values  U. 

The  solution  of  Equation  2  is  a  properly  posed  problem  if  the 
following  conditions  are  satisfied: 

(2a,)  The  solution  of  (2)  exists  for  any  u{C)  €U« 

(2b.)  The  solution  of  (2)  Is  unique  in  F(m). 

(2c.)  The  solution  of  (2)  depends  continuously  on  u(^)  in  the 
metrics  of  F  and  U. 

In  such  a  case  there  exists  a  function  Bu(C)»  defined  and  continuous 
over  all  U,  and  B  i s  an  inverse  operator  of  A,  where: 

A''u(C)  =  f(m)  =  R[m,u(C)]  (3) 

If  even  one  of  the  conditions  2a,  2b,  or  2c  is  not  satisfied  [u(^)  = 
Af(m)]  is  an  IPP.  In  such  a  case  due  to  the  unboundedness  of  the 
operator  A“1  the  solution  does  not  exist  for  all  u(^)  and  is  unstable. 
To  arbitrarily  small  variations  of  u(^)  may  correspond  arbitrarily 
large  changes  In  f(m). 

As  an  example  consider  a  Fredholm  Integral  Equation  of  the  First 
Kind. 

ttCC)  =  A  [c,f  (m)]  =  J  K  (C.  m)  f  {m)clm; 

a 


k 


c  <  Ci  d 


(4) 
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It  can  be  shown  that  the  inverse  problem  of  Equation  4,  e.g., 
computation  of  f(m)  on  the  basis  of  a  given  function  u(C  ),  is  improper. 
Indeed,  to  the  functions  f  ^  (m)  and  f2(ni)  =  f  ^  (m)  +  p  cos  aim,  f^(m), 
f2(m)  eF,  where  p  is  an  arbitrarily  large  given  number,  there  will 
correspond  functions  u^(C)  and  U2(^)  for  which  the  normllu^(^)  " 

U2(C)  II  is  arbitrari  ly  small  if  cu  is  sufficiently  large.  However, 
if  the  class  of  admissible  solutions  is  the  compact  class  F,  the  inverse 
mapping  U  ^  *F  will  be  stable  (Reference  48).  In  other  words,  for 
every  e>'0,  there  will  exist  a  S(e,  F)  such  that  the  inequality 
I  I  u^-U2  I  I  <  8  (  €  ,?)  implies  II  f  ^-f2  ll<€if  u,  and  U2  €U  =  {u(C)  = 

A[  ^  }«  where  ?”  is  a  compact  class  of  functions. 

An  IPP  will  be  called  a  Regularizable  Problem  in  Tikhonov's  sense 
if  it  has  at  least  one  regularizing  algorithm  defined  in  the  following. 
Such  regularizing  algorithms  are  practical  aids  for  solving  an  IPP. 

Let  an  IPP,  f  =  R[m,  u]  be  given,  where,  to  a  special  function  u(C), 
there  corresponds  a  solution  f (m)  =  R[m,  u(C  )]•  bet  also  an  approximating 
function  ^(C)  for  u  ( C  )  be  given,  such  that  H  u-u  ll<  8,  where  8  is 
known.  It  is  required  to  find  f (m) ,  an  approximation  to  f(m),  with  an 
assigned  precision  11  f-f  11  <  €  if  8  is  sufficiently  small. 

R[m,  u(C)]  will  be  called  a  regularizing  operator  if  the  following 
conditions  are  fulfilled: 

(1.)  Rg[m,  u(^)]  is  def  i  ned  for  every  ”u  €  U  and  S  >  0. 

(2.)  If  forlK^)  there  exists  a  f(m)  =  R[m,  u(C)],  then  for 

every  positive  e  ,  there  exists  a  S(€,f)  such  that  if 
II  iT-u  II  <  8,  then  il  f-f  II  p  <  €  where  f  g  =  R[m,  ^  . 
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Consider  a  Fredholm  Integral  Equation  of  the  First  Kind  as  expressed 

by  Equation  4,  with  the  boundaries  a=0,  b-m  ,c=0,  and  d  =  1. 

max 

If  the  kernel  K(  ^  ,m)  Is  continuous,  and  if  for  iTC^)  =  0  there  exists 
Just  one  solution  f(m)  =  0,  the  Integral  Equation  can  be  solved  In  the 
following  way  by  applying  variational  calculus  and  Tikhonov's  regularizing 
algorithm. 


Instead  of  dealing  with  the  conventional  functional: 

N[f(m)}  u  (C)]  =  /  [  A  [C,  f(m)]  -  u  (C)]  dC 

Tikhonov  suggests  applying  another  smoothing  functional: 


(5) 


Mp  [f(m)j  u(C)]  =  N[f(m),*  ^(C)]  +  aii^"^[f{m)]  (6) 

where  is  the  regularizing  functional.  Then, 

0  1  =  0 

where  the  Kj(m)'s  are  positive  continuous  functions,  and  a  is  an 
arbitrary  parameter,  which  minimizes  the  functional 


The  application  of  the  variational  Euler  Equation  and  boundary 
conditions  transform  .Equation  6  into  a  new  expression: 


i+l 

=  £  (“D  K,{m) 


i  =  0 


d'f 

dm' 


M 


•}-{  /  Klm,C)f(C)dC-b(m)}=0  (8) 
0 
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with  the  boundary  conditions 


n+l 


w  (m)  =  {  ^  (-1)'  [ K.(m)f'(nn)] 

\=i+\  ^  ^ 


(i-A-l) 


} 


=  0;  it 


0,m, 


max 


where 


_  I 

K(m,C)  =  /  K(^,m)K(C,C)df 

0 

and 

b(m)  =  J  K  (C,  m)u  (^)  dC 
0 


Equation  8  will  lead  to  an  approximate  computation  of  the 


,2, •••n  +  l)  (9) 


(lO-a) 


(lO-b) 

function  f (m) . 
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SECTION  III 
COMPUTATION 


Consider  a  Fredholm  Integral  Equation  of  the  First  Kind: 

u(C)  =  /  K(C,m) f  (m)  dm  (II) 

a 

where  K(  ^,m)  is  a  continuous  kernel  within  a  range  a<^<b  and  c<m<d. 

The  distribution  function,  f(m),  fulfills  the  following  boundary  conditions 


f(c)  =  f(d)  =  0 

and 


(12-a) 


(12-b) 


The  function  f(m)  can  be  mapped  to  generate  the  function  u(C)-  If 
an  Inverse  operation  will  be  applied  on  u(^),  a  function  f  (m)  can  be 
calculated.  Finally,  f (m)  can  be  mapped  again  to  generate  a  function 
u(C)*  The  norms  il  u(^)-u(^)  II  and  llf(m)-f(m)  II  will  show  the 
quality  of  the  transformation.  Such  a  quality  is  influenced  by  the 
following  factors: 


(a.)  The  nature  of  f{m)  function; 

(b.)  The  nature  of  the  kernel; 

(c.)  The  way  of  obtaining  u(C); 

(d.)  The  type  of  regularizing  function  applied; 

(e.)  The  type  of  matrix  used  for  approximate  calculations. 
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To  examine  these  factors,  various  distributions,  kernels,  methods  of 
integration,  regularizing  functions,  and  approximate  matrices  have  been 
cons i dered . 


The  following  f (m)  distributions  have  been  examined: 
(a-1.)  A  mono-modal  distribution: 


f  (s)  =  8s^(l-s)^,'  0<s<l  (13) 

(a-2.)  A  bi-modal  distribution: 

f  (s  )  =  f,  (s  )  +  f2(s  );  0<s<l  (14) 

where  . 

f,(s)  =  78.l25s^(s-0.6);  0<  s  <  0.6  (14-a) 

f,(s)  =  Oi  0.6  <s<l  (14-b) 

fgls)  =  Oj  0<s<0.5  {14-c) 

fgCs)  =  I62(s-0.5)*(|-s)^  0.5<s<l  (l4-d) 


The  following  kernels  have  been  applied: 

(b-1.)  The  kernel  appearing  in  Fujita's  Equation  for  concentration 
gradient,  (Equation  2)  is 

^  2 -XmC 

=  1^,-Xii.i  0«"'S"'max  i  0*^^' 

Letting  m-  X  =  0/  u(t)  =  anelT(s)  - 

f(smmax)>  Equation  15  can  be  transformed: 

- ,0<s<li  0<C<l  (15-0) 

-ps 
I  -  e 

and  Equation  2  will  be  replaced  by  Equation  16: 


u(C)  =  ' 


1  dC(£) 

""maxCo  dC 


I 

I 


K2(C,s)f  {s)ds 


(16) 
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(b-2.)  The  kernel  appearing  in  the  Fujita's  Equation  for  con¬ 
centration  (Equation  1)  is 

-XmC 

I  — '  e 

Using  a  previous  setting,  we  will  obtain: 


and 


K,(C,s)  = 


(17-0) 


C  { f ) 

^  m  C  ^  I  K,(C,8)  f  (S)d8  (18) 

max  0  Q  ® 

Two  methods  of  numerical  integration  have  been  applied: 

(c-1.)  The  integration  was  achieved  by  applying  the  trapezoidal 
rule. 

(c-2.)  The  Simpson  formula  was  applied  in  Integration. 

Various  types  of  regularizing  functions  have  been  used: 

(d-i.)n  =  iifO)  II 2 

and 

(d-2.)  a  a  =a(ll  f(’^  II  2  +  y,  II  f(2)  II  2) 


To  examine  weighting  of  elements  of  a  matrix  used  In  an  approximation 
the  variable  of  the  function  u  was  defined  in  the  two  different  ways: 


(e-l)  C, 


2  2 
^b  -f  . 
2  2  » 


and 


6-2^2  = 
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SECTION  IV 

RESULTS  OF  COMPUTATION 

To  emphasize  the  importance  of  regularization:  1)  a  mono-modal 
distribution  (Equation  13)  was  assumed,  2)  TTC^)  was  computed  by 
Equation  19  with  the  aid  of  Simpson's  formula  for  integration,  and 
3)  variational  calculus  (see  Equation  5)  without  Tikhonov's  regularizing 
function  was  applied  so  as  to  obtain  a  "best  fit"  computed  distribution. 
The  result  of  this  computation  is  shown  by  Figure  1.  Note  that  although 
the  results  are  very  erratic  the  error  analysis  II  u(^)  -u(C)  II  was 
a  minimum.  The  behavior  is  typical  for  an  IPP. 

Upon  including  a  regularizing  function  (see  Equation  7)  and 
applying  Euler's  equation  from  variational  calculus  to  the  functional 
M|[f;u]  (see  Equation  6)  a  distribution  f  “(s)  was  obtained  for  each 
value  of  the  parameter  a  .  The  objective  of  the  search  is  to  find  the 
val ue  of  a  wh i ch  wi  1  1  produce  a  minimum  error  analysis.  Using  inf  lM*^| 
as  the  error  criterion.  Figure  2  shows  a  typical  example  of  how  the 
appropriate  value  of  0!  is  determined.  This  tells  the  experimenter 
that  for  the  error  he  has  in  the  experimental  values  u(^),  this  is  the 
best  initial  approximate  distribution.  From  this  starting  approximation 
the  experimenter  can  then  apply  some  algorithm  to  infer  a  distribution 
which  will  produce  a  still  smaller  error  analysis.  He  is  now  in  a 
position  to  do  this  because  he  has  an  admissible  solution  to  his 
originally  improperly  posed  problem. 
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Using  the  concentration  kernel  (Equation  17"a)  and  the  gradient  kernel 
(Equation  15"a)  with  the  mono-modal  distribution  (Equation  13)  and  the 
regularizing  function  ^  (f) ,  the  results  obtained  are  shown  in 
Figures  3  and  k  respectively.  For  both  cases  the  Integration  involved 
In  obtaining  7(C)  was  performed  by  the  trapezoidal  rule.  Since  the 
concentration  kernel  provided  a  better  solution  to  the  problem,  especially 
for  low  molecular  weights,  this  kernel  was  used  for  all  further 
Investigations.  It  can  easily  be  shown  that  the  gradient  kernel 
(Equation  15)  degenerates  for  very  low  molecular  weights,  thus 

iifn  r  A  tn  e  i 


and  the  concentration  kernel  (Equation  17)  does  not  degenerate  under 
identical  circumstances: 

I'm  1  _  . 


Therefore,  computationally  the  concentration  kernel  Is  preferred  when 
considering  low  molecular  weights.  For  Figures  3  and  k  the  error 
criterion  minimized  was  II  u(C)"U*(C)  II  •  Here  u*(C)  was  computed 
by  using  the  calculated  function  f  (s).  To  demonstrate  the  numerical 
agreement  between  u(C)  and  u*(C)  Figure  5  presents  the  curves  for  u(C  ) 
and  [7(C)“U*(C)]  X  10^  associated  with  Figure  3*  If  and  u*(^  ) 

had  been  plotted  on  the  same  graph,  for  the  ordinate  scale  used  the 
two  curves  would  have  coalesced.  The  numerical  values  of  7(C)  and 
u*(C)  agreed  to  three  significant  figures. 


The  computation  of  a  b I -modal  distribution  proved  to  be  more  complex. 

At  first  Equation  18  was  applied  using  the  trapezoidal  rule  for  Integration. 
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The  resulting  f^(s)  obtained  by  minimizing  ||  u(^)-u*(^)  ||  was  the 
three-peaked  curve  shown  in  Figure  6.  Upon  using  Simpson's  formula  for 
integration  and  varying  a  so  as  to  minimize  II  ’u(C)“U*(C)  H  a  good 
"fit"  to  the  initial  distribution  (Equation  14)  was  obtained  (see  Figure  6). 

To  study  the  effect  of  the  error  analysis  on  choosing  the  appropriate 
value  for  CL  ,  three  criteria  were  applied: 

( I  ).  inf  I  [?(s  ),  o’  (C)]  I 

(2) .  inf  |maxlf  (s )- f°(s  )l| 

(3) .  llu(C)-u*(C)  II 

The  results  are  shown  in  Figure  7-  In  a  true  experimental  situation  the 
error  anal  ys  is  i  nf  |  max  I  f(s)-f  ^(s)  I  j-  would  not  be  ava  liable.  Choos  i  ng 
between  the  other  two  criteria  it  would  qualitatively  appear  from 
Figure  7  that  II  "u(^)-u*(^)  ||  is  better  than  one  could  expect  to 
derive  deductively  from  Tikhonov's  theory. 

For  all  computations  mentioned  above  13  =  k.Z5-  Recall  that 
/3  =  where  (3  is  proportional  to  the  square  of  the  rotor  speed 

(see  Equation  l).  If  X  =  1.02  x  and  cu  ^  =  4.1693  x  10^ 

sec  (6,166  RPM)  ,  then  X  =  4.25  x  lO”^  (mol  wt)"^.  Assuming  = 

10^,  the  value  of  f3  (“^ni^ax^  is  /3  =  4.25.  The  following  values 
for  m  were  chosen:  4.706  x  10^;  4.706  x  10^;  7*059  x  10^;  and 

m3  X 

1.8824  X  10^.  The  corresponding  values  for  ^  are  respectively  2.00, 

20.00,  30.00,  and  80.00.  The  results  obtained  by  varying  a  so  as  to 
minimize  the  error  criterion  II  u(C)-u*(C)  II  are  shown  in  Figure  8. 

No  attempt  was  made  to  determine  the  exact  optimum  value  for  (3  . 
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To  express  the  Integral  operation  (Equations  17"a  and  18)  of  the 
concentration  kernel  In  matrix  form,  the  matrix  is  weighted  by  large 
numerical  values  in  the  upper  right-hand  corner.  Since  use  of  the 
regularizing  functions  tends  to  enhance  the  diagonal  and  near-diagonal 
elements,  modification  of  the  concentration  kernel  so  as  to  have  the 
large  numerical  values  of  the  matrix  near  the  diagonal  should  greatly 
improve  the  computations.  This  was  easily  performed  by  redefining  ^ 
to  be 


Previously  C  was  defined  as 


In  using  the  same  bi-modal  distribution  applied  previously  and  the 
modified  concentration  kernel  it  was  found  that  the  computed  distribution 
f  ®(s)  was  in  poorer  agreement  with  the  assumed  f(s)  than  when  the 
unmodified  concentration  kernel  was  applied  (see  Figure  9). 
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SECTION  V 
CONCLUSIONS 

Tikhonov's  regularizing  functions  have  proved  excellent  as  a  tool  to 
solve  the  inverse  problem  associated  with  the  Fredholm  Integral  Equation 
of  the  First  Kind  which  correlates  the  ultracentrifugal  data  obtained  by 
sedimentation-diffusion  equilibrium  and  the  desired  molecular  weight 
distribution  (MWD) .  This  technique  is  especially  desirable  since  a  MWD 
can  be  obtained  from  data  at  only  one  rotor  speed.  Computationally  this 
technique  has  provided  excellent  results  when  the  MWD  is  a  mono-modal 
distribution  and  good  results  for  a  bi-modal  distribution.  However,  for 
a  multi-modal  distribution  (fine  structure)  an  additional  algorithm  has 
to  be  appl ied. 
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Figure  2 


I  I  I  >  u 

-10.2  -9.8  -9.4 


log  a 

Minimizing  the  smoothing  functional  M*^  [f  (s)  .“u  (  ^  )  ] .  A  plot 
of  log  a  versus  . 
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Figure  k.  Determination  of  MWD  for  a  mono-modai  distribution  by  applying 
Tikhonov's  regularizing  function  and  the  concentrj^t ion  kernel. 
The  solid  gradient  curve  represents  the  original  f(s)  and  the 
crosses  represent  the  computed  distribution  (s). 
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n(s:=  n/41.  5;  n=  I,  2,  — ,  41) 

Figure  5.  The  ultracentrifugal  u(^)  functions  computed  from  the 

original  distribution  and  the  calculated  distribution.  The 
solid  curve  represents  the  original  TjC^).  The  dotted  line  is 

[u(C)-uMC)]  X  103. 
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Figure  6.  Evaluation  of  the  bi-modal  distribution.  The  solid  line  represents  the  original  distribution. 
The  dotted  line  was  the  resulting  distribution  using  Simpson's  formula  for  integration. 

The  dashed  curve  resulted  from  application  of  trapezoidal  integration. 
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n(s  =  n/4l.5;  n  I,  2,  41) 

Figure  1.  Resulting  distributions  evaluated  according  to  different  error 
criterion.  The  solid  curve  is  the  original  distribution. 


AFML-TR-67-121 
Part  IV 


APPENDIX 

COMPUTER  PROGRAM  LISTING 
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IIBJOB  MAP 

JIBFTC  BEGLAR 

CGMMCN/ZYT/Xl,X2fU»A 
C  CIMENSICN  FOR  CQMMGN 

DIMENSION  U(60),A  I6C,60) 

C 

DIMENSION  XK(60«6C),S(6Q)tX(6C),ZI60)yBK(6G,60) ,B(60)tZP(60) .ALPHA 
1(10),LP(60) 


CODED  NOVEMBER  1969  BY  DoR.  WIFF  AND  M,T,  GEHATIA 
70  CHANGE  INITIAL  DISTRIBUTION,  THE  SECTION  SO  LABELED 
MUST  BE  MODIFIED, 

THIS  PROGRAM  IS  CODED  TO  RLN  ON  A  41  POINT  MESH,  IF  ANOTHER 
MESH  (EVEN  A  V/RIA8LEI  IS  tESIRED  THE  APPROPRIATE  CHANGES 
MUST  BE  MADE,  PROGRAM  SEARCHES  FOR  ALPHA  WHICH  MINUMIZES  THE 
GIVEN  ERROR  CRITERION  (SO  INDICATED),  IT  BEGINS  WITH 
loCE-LPHAl  AND  SEARCHES  TO  l,OE-LPHA2  •  THEN  THE 
MINIMUM  REGION  IS  AGAIN  SUBDIVIDED, 


REAC(5,200) 

WRITE(6,100C ) 

WRITE(6,200) 

WBITE(6,1001 ) 

CLX  =  A1.5 
DLS  =  Alo5 

XCl  =  C,78125000E  02 
XC2  =  C,162CC000E  02 
AXX  >  1, 

EVALUaTE  INITIAL  DISTRIBUTION 
DC  1  I  =  l.Al 
A1  =  I 

S(I)  »  Al/DLS 

1  X(I)  =  Al/DLX 
DC  2  I  =  1,41 
A1  =  S(n**2 

A2  =  (S( I)-l, 1**2 

2  Z(I)  =  8o*ai*A2 

XSIG  =  LAMBDA*MAXIMUM  MOLECULAR  WEIGHT 
REAC(5,103)  XSIG 
CO  4  I  =  1,41 
CCEF  =  C, 

DC  2  J  =  1,41 
A1  =  XSIG*S(J) 

A2  =  A1 

A3  =  EXP(-A1*X(  I  )  ) 

A4  =  EXP(-Al) 

A5  =  1,  -  A4 
A6  =  A2*A3/A5 
XK(I,J)  =  A6 

IF( J,EC,1,GP,J,EC«41)  GO  TO  41 
IF( JoEC.2,0R,J,EC,4C)  GO  TO  42 
KNUM  =  J/2 
JNUM  =  (J+l)/2 
IF( JNLMaEQoKNUM)  CO  TO  42 
SIG  =  2, 

GC  TC  43 

41  SIG  =  lo 
GC  TC  43 

42  SIG  =  4. 
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43  CCEF  »  CCEF  ♦  S1G*XK<  !♦  JJ’tZC  J  )/{3.*DLSI 

3  CONTIAtE 

C  CALCULATICN  OF  U(ZI»  BY  SIPPSON'S  FORMULA 

4  U«IJ  =  CCEF 
CALL  REGZCXRtBtBKtDLXI 


PRCCRAM  READS  ALPHA(NLAST»  AS  DATA  FROM  NFIRST  THROUGH 
ACCCE,  THEN  FRCGRAM  CONTINUES  FROM  NCOOE41  THROUGH 
KUPP  SEARCHING  FCR  MINIMUM  FOR  EACH  DERIVATIVE  RETAINING 
PREVIOUS  VALUES,  IF  NCCDE  =  0  ,  SEARCH  BEGINS  WITH  NFIRST, 
IF  NFLAG,GT,D  PROGRAM  READS  ONE  VALUE  OF  ALPHA  AND 
COMPUTES  FOR  ONLY  THIS  ONE  VALUE 


REAC(5,1C1J  NCODE, NFIRST, NUPP.NFLAG 
DO  30  NLAST  =  NFIRST, NUPP 
IFINFIAG,GT,C)  GO  TC  40 
IFCNLAST,LE,NCODE)  GO  TO  31 
BEACC5,100)  LPHAl,LPHA2 
LXP  =  IABS(LPHA2  -  LPHAl)  +  I 
NUM  =  C 

00  20  II  *  1,LXP 
IXP  =  LPHAl  4  II  -  1 
DO  21  KL  =  1,9 

ALPHA(NLAST)  =  FLCATIKL  ^♦10,>^■4IXP 
C  OBTAIN  HCCIFIED  MATRIX 

CALL  REG3 (NFIRST, NLAST, ALPHA, BK,B,DLX,DLS) 

C  CBTAIN  INVERSE  SOLUTION 

CALL  ECS(41) 

C  EVALUATE  ERROR 

CALL  REG4 ( Z, ZP, XK  ,DLS, UP, UAV6, NFIRST, NLAST, ALPHA, DLX I 
IF{KL,EC,1,AND,II,EC,1>  GO  TO  22 
IFrLAVG,GT,AVGll  GO  TO  28 
AVGl  =  UAVG 

C  STCRE  MINIMUM  ERROR  AND  CORRESPONDING  ALPHA 

XM  »  ALPHA(NLAST» 

NX  =  IXP 
NUM  =  C 
GO  TC  21 

22  AVGl  =  UAVG 

C  STCRE  FIRST  ALPHA  USED  AND  ASSOCIATED  ERROR 

XM  =  ALPHA(NLAST) 

NX  *  IXP 
GO  TO  21 
28  CONTINUE 
21  CCNTINLE 
20  CONTINUE 

23  CONTINUE 

XMM  =  XM  -  1C,**NX 
NUM  =  C 

DC  25  I  =  1,20 

ALPHAINLASTI  =  XMM  4  FLOAT!  I-n*10,**(NX-l» 

C  CBTAIN  MCCIFIEI  MATRIX 

CALL  REG3(NF1RST, NLAST, ALPHA, eK,B,DLX,OLS> 

C  CBTAIN  INVERSE  SOLUTION 

CALI  ECS(41I 
C  EVALUATE  ERROR 

CALL  BEG4 (Z,ZP,XK,OLS, UP, UAVG, NFIRST, NLAST, ALPHA ,DLX» 
IF(I,EC,1I  6C  TO  26 
IF(LAV6,GT,AVGl»  CO  TO  29 
AVGl  =  UAVG 
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C  STCRE  ERROR  AND  CORRESPONDING  ALPHA 

yV  =  ALPHA(MAST) 

GC  TC  25 

26  AVGl  =  UAVG 

C  STCRE  FIRST  ALPHA  USED  AND  ASSOCIATED  ERROR 

YR  =  ALPHA(ALAST) 

GC  TC  25 

29  AUM  =  NUP  +  1 
IF(MJP.EC.A)  GO  TC  27 

25  CCNTINLE 

C  START  EVALUATION  FOR  ALPHA  WITH  MINIPUH  ERROR 

27  ALPHA(NLAST)  =  YP 
GC  TC  51 

40  REAC(5,102)  ALPHA(NLAST) 

C  IF  CCPPUTATION  PROCEEDS  FOP  ONLY  ONE  ALPHA  BEGIN  HERE 

51  CCNTULE 

C  CBTAIN  MCCIFIED  PATRIX 

CALL  REG3(NFIRST,NLAST, ALPHA, EK.B.OLXtOLS) 

C  CBTAIN  INVERSE  SCLUTION 

C  EVALUATE  ERROR 

CALL  ECS(41I 

CALL  REG4( Z,ZP,XK  ,DLS , UP , UAVG  .NFIRST ,NL AST, ALPHA ,DLX ) 

DC  24  I  =  1,41 

C  Z(I)  =  ORIGINAL  CISTRieUTICN 

C  ZPd)  =  BACK  SOLUTION 

C  UCn  =  CCRRESPCNCS  TO  INPUT  DATA,  COMPUTED  USING  Z(I) 

C  LPd)  =  BACK  SOLUTION  COMPUTATION 

C  Sd)  =  VARIABLE  FOR  Zdt,  CORRESPONDING  TO 

C  PCLECULAP  WEIGHT 

WRITE (6, 2001  )  I,ZP(  II,I,Z(I),  I,UPd),I,Ud),I,Sd4 
24  CONTINUE 

WRITE(6,200C)  d,ALPHAd),I  =  NFIRST, NLAST) 

WRITE(6,2002  )  UAVG 
WRITE(6,1C4)  XSIG 
104  FCRPATdH  ,7HXSIG  =  ,1PE1C.3) 

WRITEI6,100C) 

GC  TC  30 

31  REAC(5,102)  ALPHAINLAST) 

30  CCNTINUE 
WRITE(6,1001 ) 

WRITE(6,700C  ) 

STOP 

100  FCRPaT(2I3) 

101  FGRPAT<4I2) 

102  FCRPATC 1PE9.2) 

103  FCRPAT(E10,2 ) 

200  FGRPAT(72H 

1  ) 

1000  FCRPAT(1HI/1HA) 

1001  FCRPAT (IHA/IHA ) 

2000  FCRPATdH  ,6HALPhA( ,  I2,4H)  =  ,1PE9.2) 

2001  FCRPATdH  ,7HZ-CALCI  ,  I2,4H  )  =  ,  E12.5 , 2X  ,  7HZ-TRUE  ( ,  12 ,4H )  =  ,E12a5, 
12X,7HCALC  U(,I2,4H)  *  , E12. 5 , 2X, 2HU( , 12 ,4H )  =  ,E12.5,2X,2HS( ,I2,4H 
2)  =  ,E12<.5) 

2002  FCRPATdH  ,11HM(ALPHA)  =  ,E12,5) 

7000  FCRPATdH  ,  2CX  ,  6  (  f  X  ,  lOHENC  OF  RUN)/1H1) 

END 

*IBFTC  RRR4  DECK 

SUBROUTINE  BEG4(Z  ,ZP,XK,DLS, UP, UAVG, NFIRST, NLAST, ALPHA,OLX) 

C  THIS  SUBRCUTINE  PROCESSES  THE  COMPUTED  ZPd) 

C  CALCULATES  UPd)  AND  THE  ERROR 

C  EPPCR  CRITERION 
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C0^P0A/ZYT/)Sl,X2,Ut<l 
C  DIMENSIONS  FOR  COMCN 

DIMENSION  UC60»,A(60,60) 

DIMENSION  ZP(60)»>K(60t60) fDP(60) 

DIMENSION  Z(60I,AIPHA(10) 

00  14  I  =  1,41 
ZPIll  =  All, 421 
14  CONIINUE 
UAVG  =  0. 

DC  40  I  =  1,41 
COEF  =  0, 

COEFl  =  0, 

DC  41  J  =  1,41 

IFIJ.EO.l.OR.J.EC.Al)  GO  TO  42 
IF(J.EC.2«0R.J«EC.4C)  GO  TO  43 
KNUM  =  JZ2 
JNUM  *  (J+ll/2 
lF.(JNtM,EQ,KNUMI  €0  TO  43 
SIG  =  2. 

GO  TO  44 

42  SIG  =  1, 

GO  TO  44 

43  SIG  =  4. 

44  COEF  =  COEF  +  SIG AXK 1 1 , J I ♦ ZP ( J » / (3.*DLS I 
IFCIoLT,4U  GO  TO  41 

IFIJoEC.l)  GO  TO  1 
KL  =  J-1 
GO  TO  2 

1  KL  *  1 

2  IF(J.EC.41)  GO  TO  3 
KJ  *  J*1 

GO  TO  4 

3  KJ  =  41 

4  COEFl  »  COEFl  +  S IG* ( ( ZPCK J J-ZPIKLM **2 l*DLS/3. 

41  CONTINLE 

UPCII  =  COEF 

IF( I.EC.1.0R.I.EQ.41)  GO  TO  3C 
IF(  I.EC«2,0R,UEC«40)  GO  TO  31 
KNUM  =  1/2 
JNUM  =  (I-H)/2 
IF(JNLM*EQ,KNUMJ  GO  TO  31 
SIG  =  2. 

GO  TO  23 

30  SIG  =  1, 

GO  TO  33 

31  SIG  =  4. 

33  UAVG  »  UAVG  4  SIG*( (UP C D-UC 11 )**2 I / (3. *DLX » 

40  CONTINLE 

UAV6  =  U4VG  4  ALPM(1)*C0EF1 

RETLRN 

END 

IIBFTC  RRR2  DECK 

SUBBOLTINE  RE62(XK, B, BK ,DLX ) 

THIS  SUBROUTINE  INTEGRATES  XK ( 1 , J »*XK( 1 , J J  OVER  Zl-VALUES 
TO  OBTAIN  NEW  MATRIX  6K(I,J1 
CCMMCN/ZVT/X1,X2,L,A 
DIMENSIONS  FCR  CCNMCN 
DIMENSION  U(60J,A(6C,60) 

DIMENSION  XK(60,6€),SK{60,6C),B(60) 

SIMPSON  RULE 
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CC  5  I  =  1,A1 
DC  5  J  =  ItAl 
CCEFl  *  C, 

CCEF  *  0« 

DC  2C  K  =  1,A1 

1F(K.EC.1.CF.K.EC.41)  GC  TO  21 
1F(K.EC.2.CP.K.EQ  <,4C  )  GG  TO  22 
KNUP  =  K/2 
JNUP  =  (K+l)/2 
IF( JNLP.EC.KNUM)  GO  TO  23 

22  SIG  =  2o 
GC  TC  24 

21  SIG  =  1. 

GC  TC  24 

23  SIG  =  4o 

24  A1  =  £IG^XK(K,I)*)K(K,J)/{3,*CLX» 

IFd.GT.n  GC  TC  7 

A2  =  £IG*XK(K,JI*L(K)/(3.#OLX) 

CCEFl  =  CCEFl  +  A2 

7  CCEF  =  CCEF  ♦  A1 
20  CENTIME 

IF(I.GT.l)  GC  TC  8 
eiJJ  =  CCEFl 

8  eK(I,J)  =  CCEF 
5  CCNTINLE 

RETURN 

END 

IIBFTC  RRR3  DECK 

SUBRCLTINE  R EG3 ( N FIRST , NLA  ST , ALPHA, BK,B, DUX  ,  DLS  ) 

This  SUBRCUTINE  INTRODUCES  THE  REGULARIZATION  TERMS 
IN  THE  MATRIX  EK(I,J).  TFE  FINAL  REGULARIZED  MATRIX 
IS  A(I,J| 

CCMMCN/ZYT/X1,X2, U,  A 
C  DIMENSIONS  FCR  CCFMCN 

DIMENSION  0(60) , A (60,60) 

DIMENSION  BK (60,60) ,B(60) ,ALPFA( 10) 

DC  S  I  =  1,41 
DC  S  J  =  1,41 
A(I,J)  =  BK( I,J)/GLS 

9  CCNTINLE 

DC  63  N  =  NFIRST,NLAST 

DC  60  I  =  1,41 

NUM  =  ZfN  +  1 

DC  60  K  =  1,NUM 

A2  *  CALC(N,K,DLS  ) 

A1  =  ALPFA(N)*A2 
LABEL  =  I  +  N  -  K  +  1 
64  CONTINUE 

IF(LAEELoLEoC)  GC  TC  61 
IF(LAEEL.GTo41 )  GO  TO  62 
A(I, LABEL)  =  A(I, LABEL)  +  A1 
GC  TC  60 

61  LABEL  =  lABS(LABEl)  +  1 
GO  TC  64 

62  LABEL  »  82  -  LABEL  +  1 
GC  TC  64 

60  CCNTINLE 

63  CCNTINLE 

DC  50  I  =  1,41 
SO  A(I,42)  =  B(I) 

RETURN 

END 
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$IBFTC  EEXX  DECK 

SUBPCiUKE  ECSCNI 

C  ST«KC«RD  SUBROITINE  TO  OBT4IN  THE  INVERSION  SOLUTION,  THE 

C  SOLUTION  IS  IK  COLUMN  J  +  l  IN  Ad.JI  MATRIX, 

CCMMCK/ZVTXX1,X2,U,A 
C  CIMENSICN  FOB  COMMON 

DIMENSION  UCEOl.A  (6C,60I 
DIMENSION  BI60) 

53  Nl=N4l 
5  M»N1 

13  00  23  K=1,N 

Kl*K4l 

lF(A(K,Kn2Ctl4,2C 

14  IF(K-K 115,19,15 

15  DC  18  l=Kl,N 
IFIA(L,Km6,18,H 

16  DC  IT  J=K,M 
B(JJ=A IK,J» 

A(K,J)  =  Aa,J) 

17  AIL,J)=-e(J> 

GO  TO  20 

18  CCNTIM'E 

19  WR1TE(6,1C1) 

GO  TC  30 

20  DC  23  1=1, N 
lFn-KI2l,23,21 

21  C=A(1,K)/A(K,K) 

DC  22  J=1,M 

22  A(I,J)=A(I,J)  -  C»A(K,J) 

23  CCNTINLE 

46  DO  29  J=M,M 
DC  25  1=1, N 

25  A(I,J)  =  A(1,JI/A(1  ,n 
29  CCNTINLE 

101  FCBMAT(19H  MATRIX  IS  SINGULAR! 

30  RETURN 
END 

tlBFTC  CALF  DECK 

FUNCTION  CALC(N,K,OLS» 

C  THIS  FUNCTION  SUBROUTINE  EVALUATES  THE 

C  CCEFFICIENTS  (BINOMIAL),  ETC, 

C  CALLED  BV  REG3 

C 

L  =  24N 
M  =  L  -  K  4  1 
J  «  K  -  1 

IF(K,EC,l,CR,K.EC,L4l)  GO  TO  10 


11 

s 

1 

12 

= 

1 

13 

s 

1 

DC 

1 

I  =  l,L 

1 

11 

s 

11*1 

00 

2 

1  =  1,M 

2 

12 

= 

12*1 

DC 

3 

I  =  1,J 

3 

13 

3 

13*1 

14 

= 

11/112*13) 

15 

3 

(-1)**(N4K4I ) 

X  =  FLCATIH’AIS) 
CALC  =  X#CLS**L 
GC  TC  99 
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10  lA  =  (-1)**N 
A1  =  FLOATdA) 

CALC  =  A1=»CLS**L 
99  RETURN 
END 

$CATA 

CONCENTRATICN  KERNEL. < IFTS. , 8# (V* (V- 1 ) >**2 
C.425E  01 
0  111 
1,00E-C7 
SEOF 
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